#!/bin/tcsh -ef
####
#
#############################################################################
#                                                                           #
# Title: Sharpen Refined Map                                                #
#                                                                           #
# (C) 2dx.org, GNU Plublic License.                                         #
#                                                                           #
# Created..........: 21/01/2015                                             #
# Last Modification: 21/01/2015                                             #
# Author...........: Nikhil Biyani                                          #
#                                                                           #
#############################################################################
#
# SORTORDER: 90
#
# MANUAL: This script takes the refined map and sharpens it using the reference's density histogram and structure factors
#
#
# PUBLICATION: 3D reconstruction of two-dimensional crystals: <A HREF="http://www.ncbi.nlm.nih.gov/pubmed/26093179">Arch Biochem Biophys 581, 68-77 (2015)</A>
# PUBLICATION: 3D Reconstruction from 2D crystal image and diffraction data: <A HREF="http://dx.doi.org/10.1016/S0076-6879(10)82004-X">Methods Enzymol. 482, Chapter 4, 101-129 (2010)</A>
# PUBLICATION: 2dx - Automated 3D structure reconstruction from 2D crystal data: <A HREF="http://journals.cambridge.org/action/displayAbstract?aid=1943200">Microscopy and Microanalysis 14(Suppl. 2), 1290-1291 (2008)</A>
# PUBLICATION: 2dx_merge - Data management and merging for 2D crystal images: <A HREF="http://dx.doi.org/10.1016/j.jsb.2007.09.011">J. Struct. Biol. 160(3), 375-384 (2007)</A>
#
# DISPLAY: SYM
# DISPLAY: realcell
# DISPLAY: realang
# DISPLAY: ALAT
# DISPLAY: number_of_beads
# DISPLAY: density_threshold_bead
# DISPLAY: maximum_resolution_reference_sharpening
# DISPLAY: sharpening_algo
# DISPLAY: number_sharpening_iterations
# DISPLAY: sharpening_reference_PDB_file
# DISPLAY: RESMIN
# DISPLAY: RESMAX
# DISPLAY: calculate_subvolume
# DISPLAY: bfactor_sharpening
# DISPLAY: if_bfactor_sharpen
#
#$end_local_vars
#
set bin_2dx = ""
set proc_2dx = ""
#
set SYM = ""
set realcell = ""
set realang = ""
set ALAT = ""
set RESMIN = ""
set RESMAX = ""
set calculate_subvolume = ""
set number_of_beads = ""
set density_threshold_bead = ""
set maximum_resolution_reference_sharpening = ""
set sharpening_algo = ""
set number_sharpening_iterations = ""
set sharpening_reference_PDB_file = ""
set bfactor_sharpening = ""
set if_bfactor_sharpen = ""
#
#$end_vars
#
set scriptname = 2dx_sharpenRefinedMap
#
set split = ($realcell:as/,/ /)
set cellx = $split[1]
set celly = $split[2]
#
echo "cellx = ${cellx}"
echo "celly = ${celly}"
echo "cellz = ${ALAT}"
#
set date = `date`
echo date = ${date}
#
\rm -f LOGS/${scriptname}.results
#
set ccp4_setup = 'y'
source ${proc_2dx}/initialize
#
#
set refined_map = "refined.map"
#
if (! -f ${refined_map} ) then
  ${proc_2dx}/protest "Refined map not found! First, do a round of RefineMergedMap"
endif
#
#
#
echo "<<@progress: 5>"
#
if(${sharpening_algo} == "1") then
    ###########################################################################
    ${proc_2dx}/linblock "Preparing appropriate files for bead-model map"
    ###########################################################################
    #
    set bead_model_map = "bead_model.map"
    #
    \rm -f ${bead_model_map}
    #---------------------------------------------------------------------------
    echo ":Launching ${bin_2dx}/create_bead_model.exe --mrcin ${refined_map} --mrcout ${bead_model_map} -b ${number_of_beads} --threshold ${density_threshold_bead} --res ${maximum_resolution_reference_sharpening}"
    #---------------------------------------------------------------------------
    ${bin_2dx}/create_bead_model.exe --mrcin ${refined_map} --mrcout ${bead_model_map} -b ${number_of_beads} --threshold ${density_threshold_bead} --res ${maximum_resolution_reference_sharpening}
    #
    echo "# IMAGE: ${bead_model_map} <Bead model map>" >> LOGS/${scriptname}.results
    #
    set bead_model_extended_map = "bead_model_extended.map"
    set bead_model_sub_map = "bead_model_sub.map"
    #
    \rm -f ${bead_model_extended_map}
    #
    source ${proc_2dx}/2dx_extend_map.com ${bead_model_map} ${bead_model_extended_map}
    #
    echo "# IMAGE-IMPORTANT: ${bead_model_extended_map} <Bead model extended map 2X2X1 unit cells>" >> LOGS/${scriptname}.results
    #
    \rm -f ${bead_model_sub_map}
    #
    if ( ${calculate_subvolume}x != "0x" ) then 
        source ${proc_2dx}/2dx_create_subvolume.com ${calculate_subvolume} ${bead_model_extended_map} ${realcell} ${ALAT} ${bead_model_sub_map}
        #
        echo "# IMAGE-IMPORTANT: ${bead_model_sub_map} <Bead model sub map>" >> LOGS/${scriptname}.results
    endif
    #
    set reference_map = ${bead_model_map}
endif
#
if(${sharpening_algo} == "2") then
    #------------------------------------------------------------------------
    ${proc_2dx}/linblock "E2PDB2MRC.PY: Converting pdb to mrc"
    #------------------------------------------------------------------------
    #
    set model_e2pdb2mrc = "SCRATCH/model_e2pdb2mrc.mrc"
    #
    \rm -f ${model_e2pdb2mrc}
    #
    e2pdb2mrc.py -R ${maximum_resolution_reference_sharpening} -A 1.0 -B ${ALAT} ${sharpening_reference_PDB_file} ${model_e2pdb2mrc}
    #
    #------------------------------------------------------------------------
    ${proc_2dx}/linblock "E2PROC3D.PY: Correcting the size of bead model mrc"
    #------------------------------------------------------------------------
    #
    set model_mrc = "model.mrc"
    #
    \rm -f ${model_mrc}
    #
    set size_x = `printf %.0f ${cellx}`
    set size_y = `printf %.0f ${celly}`
    set size_z = `printf %.0f ${ALAT}`
    set input_size = "${size_x},${size_y},${size_z}"
    echo "Chopping to ${input_size}"
    #
    e2proc3d.py --clip ${input_size} ${model_e2pdb2mrc} ${model_mrc}
    #
    set model_map = "model.map"
    \rm -f ${model_map}
    #
    ${bin_2dx}/2dx_processor.exe --mrcin ${model_mrc} --mrcout ${model_map}
    #
    echo "# IMAGE-IMPORTANT: ${model_map} <MAP: PDB model density map>" >> LOGS/${scriptname}.results
    #
    set reference_map = ${model_map}
    #
endif
#
echo "<<@progress: +30>"
#
###########################################################################
${proc_2dx}/linblock "Preparing appropriate files for sharpening the map"
###########################################################################
#
set sharpened_hkl = "sharpened_LeftHanded.hkl"
set sharpened_map = "sharpened.map"
\rm -f ${sharpened_hkl}
\rm -f ${sharpened_map}
#
touch SCRATCH/sharpened_dummy
\rm -f SCRATCH/sharpened_*
#
#
if(${sharpening_algo} != "0") then
    #------------------------------------------------------------------------------
    echo ":Launching ${bin_2dx}/2dx_sf_hist_refiner.exe --mrcin ${refined_map} --refin ${reference_map} --temp SCRATCH/ --res ${RESMAX} -s ${SYM} --iterations ${number_sharpening_iterations} --hklout ${sharpened_hkl} --mrcout ${sharpened_map}"
    #------------------------------------------------------------------------------
    #
    ${bin_2dx}/2dx_sf_hist_refiner.exe --mrcin ${refined_map} --refin ${reference_map} --temp SCRATCH/ --res ${RESMAX} -s ${SYM} --iterations ${number_sharpening_iterations} --hklout ${sharpened_hkl} --mrcout ${sharpened_map}
    #
    set num = 1
    while ( ${num} <= ${number_sharpening_iterations} ) 
      echo "# IMAGE: SCRATCH/sharpened_${num}.map <MAP: Sharpened scratch map, iteration ${num}>" >> LOGS/${scriptname}.results
      @ num += 1
    end
    #------------------------------------------------------------------------------
    #echo ":Launching ${bin_2dx}/2dx_processor.exe --mrcin ${sharpened_map} --mrcout ${sharpened_map} --res ${RESMAX}"
    #------------------------------------------------------------------------------
    #${bin_2dx}/2dx_processor.exe --mrcin ${sharpened_map} --mrcout ${sharpened_map} --res ${RESMAX}
    #
else if(${if_bfactor_sharpen} == "y") then
    cp ${refined_map} ${sharpened_map}
endif
#
#
if(${if_bfactor_sharpen} == "y") then
    #------------------------------------------------------------------------------
    echo ":Launching ${bin_2dx}/2dx_processor.exe --mrcin ${sharpened_map} --hklout ${sharpened_hkl} --mrcout ${sharpened_map} --res ${RESMAX} --bfactor ${bfactor_sharpening}"
    #------------------------------------------------------------------------------
    ${bin_2dx}/2dx_processor.exe --mrcin ${sharpened_map} --hklout ${sharpened_hkl} --mrcout ${sharpened_map} --res ${RESMAX} --bfactor ${bfactor_sharpening}
    #
endif
#
#
if(${sharpening_algo} != "0" || ${if_bfactor_sharpen} == "y") then
    echo "# IMAGE: ${sharpened_hkl} <Sharpened HKL (MRC lefthanded) [H K L AMP PHASE FOM]>" >> LOGS/${scriptname}.results
    echo "# IMAGE: ${sharpened_map} <Sharpened map>" >> LOGS/${scriptname}.results
    #
    #
    set sharpened_mtz = "merge3Dref_Sharpened_MRClefthanded.mtz"
    \rm -f ${sharpened_mtz}
    source ${proc_2dx}/2dx_hkl_to_mtz.com ${sharpened_hkl} ${realcell} ${ALAT} ${realang} ${RESMIN} ${RESMAX} ${sharpened_mtz}
    echo "# IMAGE-IMPORTANT: ${sharpened_mtz} <MTZ: Sharpened Reference 3D MTZ file (MRC lefthanded) [H,K,L,F,P,FOM,SIGF] >" >> LOGS/${scriptname}.results
    #
    set sharpened_extended_map = "sharpened_extended.map"
    set sharpened_sub_map = "sharpened_sub.map"
    #
    \rm -f ${sharpened_extended_map}
    #
    source ${proc_2dx}/2dx_extend_map.com ${sharpened_map} ${sharpened_extended_map}
    #
    echo "# IMAGE-IMPORTANT: ${sharpened_extended_map} <Sharpened extended map 2X2X1 unit cells>" >> LOGS/${scriptname}.results
    #
    \rm -f ${sharpened_sub_map}
    #
    if ( ${calculate_subvolume}x != "0x" ) then 
        source ${proc_2dx}/2dx_create_subvolume.com ${calculate_subvolume} ${sharpened_extended_map} ${realcell} ${ALAT} ${sharpened_sub_map}
        #
        echo "# IMAGE-IMPORTANT: ${sharpened_sub_map} <Sharpened sub map>" >> LOGS/${scriptname}.results
    endif
    #
endif
echo "<<@progress: +20>"
#
#
#############################################################################
#############################################################################
${proc_2dx}/linblock "${scriptname} normal end."
#############################################################################
#
echo "<<@progress: 100>>"
#
exit
#
# These are here only to make sure they show up in the GUI:
source ${proc_2dx}/2dx_hkl_to_mtz.com
source ${proc_2dx}/2dx_extend_map.com
source ${proc_2dx}/2dx_create_subvolume.com
#
